EDA_NewZealand

Get map and extract crs for the map

nz.map <- ne_countries(country = 'New Zealand', returnclass = "sf", scale = 'medium')
nz.crs <- crs(nz.map)

nz_cropped <- st_crop(nz.map, xmin = 130, xmax = 180,
                                    ymin = -70, ymax = -33)

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

ggplot() + geom_sf(data = nz_cropped) + theme_bw()

# Get the coastline from the rnaturalearth package for comparison
nz <- ne_countries(scale = "medium", returnclass = "sf", country = "New Zealand") %>%
  st_transform(crs = proj_nzsf()) %>%
  st_crop(get_statistical_areas(area = "EEZ"))

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

Load and pre-process catalogue

Let’s look at the raw catalogue

head ../../../data/NewZealand/hypnzw23R2001_2011.xyze

##   YR MODY HRMN  SEC  LATITUDE LONGITUDE  DEPTH   MAG   NO RMSRES     x        y    GAP DMIN RZDM NP  NS SE_OT  SE_H  SE_Z  Q
## 2001 0101 0048 35.29 -39.5230  175.6982  -0.42  2.00    6  0.22    26.16  -342.11  351  33  0.0   3   3  0.24 16.17 23.37  E
## 2001 0101 0153 46.26 -42.2526  172.6404   7.57  2.60   10  0.07    18.33    55.54  172  59  0.1   5   5  0.06  1.00  2.00  C
## 2001 0101 0224 43.34 -38.0645  178.6776  35.40  2.80    6  0.02   133.69  -628.27  290  38  0.9   3   3  0.07  2.00  4.00  E
## 2001 0101 0301 37.34 -38.0325  178.7334  29.31  3.10    7  0.05   135.52  -634.06  255  42  0.7   4   3  0.14  1.50  3.00  D
## 2001 0101 0311 17.49 -38.0626  178.7789  30.95  2.90    7  0.07   140.73  -633.74  265  46  0.7   4   3  0.10  1.50  3.00  D
## 2001 0101 0318 23.15 -41.2742  172.6529 151.57  3.50   22  0.11   -51.04   -28.06  112  53  2.9  16   6  0.08  0.60  0.85  A
## 2001 0101 0403 21.72 -39.2606  175.4920   6.87  1.30    8  0.11    -5.86  -353.93  302   5  1.4   4   4  0.12  2.00  4.00  E
## 2001 0101 0437  9.63 -40.2849  173.2374 187.26  3.20    4  0.52   -83.92  -143.86  253  63  3.0   2   2  0.10 10.79 20.00  E
## 2001 0101 0637 19.56 -39.2356  175.5181  11.32  2.00    8  0.09    -5.82  -357.51  186   4  2.8   7   1  0.16  1.00  2.00  C

To parse the datetimes, I first pre-process the catalogue to split the MODY and HRMN columns. This makes it easier to parse later on dealine with these as integers etc. Alternatively, I could have loaded as strings.

Now we can load the catalogue in R and parse the date-times

# Load the reformatted catalogue
df_cat <- read.table("../../../data/NewZealand/hypnzw23R2001_2011_MNreformat.xyze", header=TRUE, sep="", col.names=c("YR","MO","DY","HR","MN","SEC","LATITUDE","LONGITUDE","DEPTH","MAG","NO","RMSRES","x","y","GAP","DMIN","RZDM","NP","NS","SE_OT","SE_H","SE_Z","Q"), fill=TRUE)

# add event number
df_cat$event_num <- seq.int(nrow(df_cat))

# parse date-times and set appropriate digital second resoluton
op <- options(digits.secs=2)
options(op)
df_cat$date <- ymd( paste(df_cat$YR, df_cat$MO, df_cat$DY, sep="-") )
df_cat$dateTime <- ymd_hms( paste(df_cat$YR,"-", df_cat$MO,"-", df_cat$DY," " , df_cat$HR,"-",df_cat$MN,"-" ,df_cat$SEC, sep="") )

## Warning: 11 failed to parse.

# sp version of the catalogue
df_cat.sp <- data.frame(df_cat)
sp::coordinates(df_cat.sp) <- c("LONGITUDE", "LATITUDE")
crs_wgs84 <- CRS(SRS_string='EPSG:4326')
proj4string(df_cat.sp) <- crs_wgs84

# sf version of the catalogue
df_cat.sf<- st_as_sf(df_cat.sp)

head(df_cat)

##     YR MO DY HR MN   SEC LATITUDE LONGITUDE  DEPTH MAG NO RMSRES      x       y
## 1 2001  1  1  0 48 35.29 -39.5230  175.6982  -0.42 2.0  6   0.22  26.16 -342.11
## 2 2001  1  1  1 53 46.26 -42.2526  172.6404   7.57 2.6 10   0.07  18.33   55.54
## 3 2001  1  1  2 24 43.34 -38.0645  178.6776  35.40 2.8  6   0.02 133.69 -628.27
## 4 2001  1  1  3  1 37.34 -38.0325  178.7334  29.31 3.1  7   0.05 135.52 -634.06
## 5 2001  1  1  3 11 17.49 -38.0626  178.7789  30.95 2.9  7   0.07 140.73 -633.74
## 6 2001  1  1  3 18 23.15 -41.2742  172.6529 151.57 3.5 22   0.11 -51.04  -28.06
##   GAP DMIN RZDM NP NS SE_OT  SE_H  SE_Z Q event_num       date
## 1 351   33  0.0  3  3  0.24 16.17 23.37 E         1 2001-01-01
## 2 172   59  0.1  5  5  0.06  1.00  2.00 C         2 2001-01-01
## 3 290   38  0.9  3  3  0.07  2.00  4.00 E         3 2001-01-01
## 4 255   42  0.7  4  3  0.14  1.50  3.00 D         4 2001-01-01
## 5 265   46  0.7  4  3  0.10  1.50  3.00 D         5 2001-01-01
## 6 112   53  2.9 16  6  0.08  0.60  0.85 A         6 2001-01-01
##              dateTime
## 1 2001-01-01 00:48:35
## 2 2001-01-01 01:53:46
## 3 2001-01-01 02:24:43
## 4 2001-01-01 03:01:37
## 5 2001-01-01 03:11:17
## 6 2001-01-01 03:18:23

EDA for the whole New Zealand catalogue

Histogram of event depths

ggplot(df_cat, aes( x=DEPTH )) + 
  geom_histogram( binwidth=1 ) + 
  coord_flip (xlim = c(350, 0) ) +
  ggtitle("Hight res. histogram of event depths for depth inversion aretfacts")

## Plot of event magnitudes using natural time

ggplot(df_cat, aes(x=event_num, y=MAG))  + geom_point(size = 0.1) 

Plot of event magnitudes using date-time

ggplot(df_cat, aes(x=dateTime, y=MAG)) + 
  geom_point(size = 0.1) +
  ggtitle("New Zealand magnitude time series")

## Warning: Removed 11 rows containing missing values (geom_point).

# Filtered for M>4
ggplot(df_cat[df_cat$MAG>4,], aes(x=dateTime, y=MAG)) + 
  geom_point(size = 0.1) +
  ggtitle("New Zealand magnitude timeseries for M>4")

Frequency-magnitude analysis for whole New Zealand

minMag <- 1
maxMag <- max(df_cat$MAG)

mags <- df_cat[df_cat$MAG>minMag,]$MAG

tmp <- hist(mags, breaks=seq(minMag-0.05,maxMag+0.1,0.1), plot=FALSE)

N.counts <- length( tmp$counts)
tmp$cumulativeCounts <- cumsum(tmp$counts[N.counts:1])[N.counts:1]

m.min <- 4
bin_m.min <- which(tmp$mids==m.min)
freq_m.min <- tmp$counts[bin_m.min]
b <- 1.1
x <- tmp$mids
y <- freq_m.min * 10^(-b*(x-m.min))
y.cum <- tmp$cumulativeCounts[bin_m.min] * 10^(-b*(x-m.min))

ggplot() +
  geom_point( aes(x=tmp$mids, y=tmp$counts) ) +
  geom_point( aes(x=tmp$mids, y=tmp$cumulativeCounts) , color='red', pch="+", size=2) +
  scale_y_log10() +
  ggtitle(paste("Frequency-magnitude plot with arbitary GR dist: b =", b)) +
  xlab("Magnitude") +
  ylab("log10(Frequency)") +
  geom_line(aes(x=x, y=y)) +
  geom_line(aes(x=x, y=y.cum), color='red') +
  geom_vline( xintercept=m.min, lty=2 )

## Warning: Transformation introduced infinite values in continuous y-axis

b.stability.list <- c()
b.error.list <- c()
m.mean <- c()

max.index.x <- length(x)-5

for( i in 1:max.index.x ){
  mag.threshold <- x[i]
  N.events <- length(mags[mags > mag.threshold])
  m.mean <- mean( mags[mags > mag.threshold] )
  m.sd <- sd(mags[mags > mag.threshold])
  b.stability.list[i] <- log10( exp(1) ) / (m.mean - mag.threshold - 0.05)
  b.error.list[i] <- 2.3 * b.stability.list[i]^2 * m.sd / (sqrt(N.events*(N.events-1)))
}

ggplot() +
  geom_line( aes(x=x[1:max.index.x], y=b.stability.list) ) +
  geom_line( aes(x=x[1:max.index.x], y=b.stability.list+b.error.list), color=2, lty=2 ) +
  geom_line( aes(x=x[1:max.index.x], y=b.stability.list-b.error.list), color=2, lty=2 ) +
  xlab("Magnitude threshold") +
  ylab("Aki b-value estimate") +
  geom_hline(yintercept = 1, lty=3) +
  ggtitle("b-value stability plot for New Zealand")

Various maps of the events

ggplot() +
  geom_hex(data = df_cat[df_cat$MAG>4,], aes(x = LONGITUDE, y = LATITUDE), bins = 50) +
  scale_fill_continuous(type = "viridis") +
  geom_sf(data = nz_cropped, fill=alpha("lightgrey", 0), color = 'orange', size=0.2) + 
  ggtitle("Density plot for M>4 events") +
  theme_bw()

ggplot() +
  geom_sf(data = df_cat.sf[df_cat$MAG>4,], size = 0.05) +
  geom_sf(data = nz_cropped, fill=alpha("lightgrey", 0), color = 'green') +
  geom_sf(data = df_cat.sf[df_cat$MAG>6,], size = 0.5, color='orange') +
  geom_sf(data = df_cat.sf[df_cat$MAG>7,], size = 0.5, color='red') +
  ggtitle("Map of event locations")

ggplot() +
  geom_sf(data = nz, fill = alpha("lightgrey", 0), colour = "green") +
  geom_sf(data = df_cat.sf[df_cat$MAG>4,], size = 0.05) +
  geom_sf(data = nz, fill = alpha("lightgrey", 0), colour = "green") +
  annotation_north_arrow(location = "tl", which_north = "grid", style = north_arrow_orienteering, height=unit(1., "cm"), width=unit(1., "cm")) +
  annotation_scale(location = "br", unit_category = "metric") +
  geom_sf(data = df_cat.sf[df_cat$MAG>6,], size = 0.5, color='orange') +
  geom_sf(data = df_cat.sf[df_cat$MAG>7,], size = 0.5, color='red', show.legend=points) +
  ggtitle("Map of event locations")

Latitude-time plot for whole New Zealand

ggplot() + 
  geom_point(data=df_cat[df_cat$MAG>4,], aes(dateTime, LATITUDE), size=0.1) +
  geom_point(data=df_cat[df_cat$MAG>6,], aes(dateTime, LATITUDE), size=1.2, color='orange') +
  geom_point(data=df_cat[df_cat$MAG>7,], aes(dateTime, LATITUDE), size=1.5, color='red') +
  ggtitle("New Zealand latitude-time plot")

Extract the South Island

#eventDate <- ymd_hms( '2009-07-15 09:22:29' )

latLims <- c( -42, -47)
minMAG <- 4

# Subset the main catalogue
df_cat.south <- df_cat[df_cat$MAG >= minMAG, ]
df_cat.south <- df_cat.south[ (df_cat.south$LATITUDE<latLims[1]), ]
df_cat.south <- df_cat.south[ (df_cat.south$LATITUDE>latLims[2]), ]
#df_cat.subset <- df_cat.subset[ (df_cat.subset$LONGITUDE>longLims[1]), ]
#df_cat.subset <- df_cat.subset[ (df_cat.subset$LONGITUDE<longLims[2]), ]

endDate <- max(df_cat.south$dateTime)
startDate <- min(df_cat.south$dateTime)
longLims <- c( min(df_cat.south$LONGITUDE), max(df_cat.south$LONGITUDE) )

head(df_cat.south)

##        YR MO DY HR MN   SEC LATITUDE LONGITUDE  DEPTH MAG NO RMSRES      x
## 163  2001  1  4 20  9 58.04 -45.0396  167.3916  91.06 4.0 18   0.14 -88.59
## 796  2001  1 19  6 54 23.42 -43.6561  169.6835   5.00 4.0 22   0.14 -60.24
## 1067 2001  1 22 20 41 24.25 -44.9499  167.4568 101.89 4.6 21   0.09 -91.82
## 1496 2001  1 31  0 57 22.26 -43.4153  179.7593  36.00 4.0 43   0.47 557.52
## 1552 2001  2  1  5  2 29.32 -44.5924  168.1326  76.88 4.0 21   0.13 -80.52
## 2190 2001  2 12  3  6 43.91 -46.7453  165.1765  12.75 4.3  7   0.35 -76.92
##            y GAP DMIN RZDM NP NS SE_OT  SE_H  SE_Z Q event_num       date
## 163   569.18 187   52  1.8 12  6  0.08  1.00  2.00 C       163 2001-01-04
## 796   331.81 139   35  0.1 17  5  0.02  1.00  5.00 C       796 2001-01-19
## 1067  558.42 191   65  1.6 13  8  0.06  1.00  2.00 C      1067 2001-01-22
## 1496 -199.50 305  433  0.1 39  4  0.32  2.38  4.00 E      1496 2001-01-31
## 1552  492.64 209   86  0.9 14  7  0.05  1.00  2.00 C      1552 2001-02-01
## 2190  825.47 315  244  0.1  6  1  1.03 21.36 26.15 E      2190 2001-02-12
##                 dateTime
## 163  2001-01-04 20:09:58
## 796  2001-01-19 06:54:23
## 1067 2001-01-22 20:41:24
## 1496 2001-01-31 00:57:22
## 1552 2001-02-01 05:02:29
## 2190 2001-02-12 03:06:43

ggplot() +
  geom_sf(data = df_cat.sf[df_cat$MAG>4,], size = 0.05) +
  geom_sf(data = nz_cropped, fill=alpha("lightgrey", 0), color = 'green') +
  geom_sf(data = df_cat.sf[df_cat$MAG>6,], size = 0.5, color='orange') +
  geom_sf(data = df_cat.sf[df_cat$MAG>7,], size = 0.5, color='red') +
  ggtitle("Map of event locations")+
  geom_rect(aes(xmin=longLims[1], xmax=longLims[2], ymin=latLims[1], ymax=latLims[2]), color="blue",fill=NA)

ggplot() + 
  geom_point(data=df_cat[df_cat$MAG>4,], aes(dateTime, LATITUDE), size=0.1) +
  geom_point(data=df_cat[df_cat$MAG>6,], aes(dateTime, LATITUDE), size=1.2, color='orange') +
  geom_point(data=df_cat[df_cat$MAG>7,], aes(dateTime, LATITUDE), size=1.5, color='red') +
  ggtitle("New Zealand latitude-time plot with Southern Island highlighted") +
  geom_rect(aes(xmin = as.POSIXct(ymd("2001-1-1")), xmax = as.POSIXct(ymd("2012-1-1")), ymin = latLims[1], ymax = latLims[2]), alpha = 0.4, fill='blue')

ggplot(df_cat.south, aes(x=dateTime, y=MAG)) + 
  geom_point(size = 0.1) +
  ggtitle("Southern Island magnitude timeseries for M>4") +
  geom_rect( aes(xmin = as.POSIXct(startDate), xmax = as.POSIXct(endDate), ymin = minMAG, ymax = max(df_cat.south$MAG+0.2)), alpha = 0.4, fill=NA, color="blue" )

## Warning: Use of `df_cat.south$MAG` is discouraged. Use `MAG` instead.

minMag <- 4
maxMag <- max(df_cat.south$MAG)

mags <- df_cat.south[df_cat.south$MAG>=minMag,]$MAG

tmp <- hist(mags, breaks=seq(minMag-0.05,maxMag+0.1,0.1), plot=FALSE)

N.counts <- length( tmp$counts)
tmp$cumulativeCounts <- cumsum(tmp$counts[N.counts:1])[N.counts:1]

m.min <- 4
bin_m.min <- which(tmp$mids==m.min)
freq_m.min <- tmp$counts[bin_m.min]
b <- 1.1
x <- tmp$mids
y <- freq_m.min * 10^(-b*(x-m.min))
y.cum <- tmp$cumulativeCounts[bin_m.min] * 10^(-b*(x-m.min))

ggplot() +
  geom_point( aes(x=tmp$mids, y=tmp$counts) ) +
  geom_point( aes(x=tmp$mids, y=tmp$cumulativeCounts) , color='red', pch="+") +
  scale_y_log10() +
  ggtitle(paste("Frequency-magnitude plot with arbitary GR dist: b =", b)) +
  xlab("Magnitude") +
  ylab("log10(Frequency)") +
  geom_line(aes(x=x, y=y)) +
  geom_line(aes(x=x, y=y.cum), color='red') +
  geom_vline( xintercept=m.min, lty=2 )

## Warning: Transformation introduced infinite values in continuous y-axis

b.stability.list <- c()
b.error.list <- c()
m.mean <- c()

max.index.x <- length(x)-5

for( i in 1:max.index.x ){
  mag.threshold <- x[i]
  m.mean[i] <- mean( mags[mags > mag.threshold], na.rm=TRUE )
  b.stability.list[i] <- log10( exp(1) ) / (m.mean[i] - mag.threshold - 0.05)
  b.error.list[i] <- 2.3 * b.stability.list[i]^2 * sd(mags[mags > mag.threshold]) / (sqrt(N.events*(N.events-1)))
}

ggplot() +
  geom_line( aes(x=x[1:max.index.x], y=b.stability.list) ) +
  geom_line( aes(x=x[1:max.index.x], y=b.stability.list+b.error.list), color=2, lty=2 ) +
  geom_line( aes(x=x[1:max.index.x], y=b.stability.list-b.error.list), color=2, lty=2 ) +
  xlab("Magnitude threshold") +
  ylab("Aki b-value estimate") +
  geom_hline(yintercept = 1, lty=3) +
  ggtitle("b-value stability plot for Southern Island")

Extract the Northern Island

#eventDate <- ymd_hms( '2009-07-15 09:22:29' )

latLims <- c( -34, -41)
minMAG <- 4

# Subset the main catalogue
df_cat.north <- df_cat[df_cat$MAG >= minMAG, ]
df_cat.north <- df_cat.north[ (df_cat.north$LATITUDE>latLims[2]), ]
#df_cat.north <- df_cat.north[ (df_cat.north$LATITUDE<latLims[2]), ]
#df_cat.subset <- df_cat.subset[ (df_cat.subset$LONGITUDE>longLims[1]), ]
#df_cat.subset <- df_cat.subset[ (df_cat.subset$LONGITUDE<longLims[2]), ]

endDate <- max(df_cat.north$dateTime)
startDate <- min(df_cat.north$dateTime)
longLims <- c( min(df_cat.north$LONGITUDE), max(df_cat.north$LONGITUDE) )

head(df_cat.north)

##       YR MO DY HR MN   SEC LATITUDE LONGITUDE  DEPTH MAG NO RMSRES      x
## 14  2001  1  1 11 58 50.85 -38.0754  178.7195  27.32 4.0 25   0.08 137.38
## 39  2001  1  2  5 34 42.69 -39.0541  176.8003  38.82 4.5 46   0.29  69.16
## 42  2001  1  2  6 22 16.54 -40.5215  175.6788  37.57 4.4 51   0.14  93.54
## 150 2001  1  4 15 46 53.47 -38.0652  178.7116  33.04 4.8 46   0.15 136.14
## 180 2001  1  5  5 30 45.62 -38.4291  175.6323 179.26 4.1 31   0.20 -53.64
## 203 2001  1  5 17  4 38.86 -38.6254  175.1182 234.35 4.1 39   0.15 -75.15
##           y GAP DMIN RZDM NP NS SE_OT SE_H SE_Z Q event_num       date
## 14  -629.49 252   31  0.9 18  7  0.06 1.50 3.00 D        14 2001-01-01
## 39  -441.69  82   49  0.8 41  5  0.08 0.41 1.79 B        39 2001-01-02
## 42  -254.02  74   18  2.1 40 11  0.02 0.24 0.32 A        42 2001-01-02
## 150 -629.99 249   40  0.8 36 10  0.08 1.50 3.00 D       150 2001-01-04
## 180 -433.87 186   64  2.8 24  7  0.12 1.00 2.00 C       180 2001-01-05
## 203 -388.89 144   35  6.7 31  8  0.12 0.66 1.03 B       203 2001-01-05
##                dateTime
## 14  2001-01-01 11:58:50
## 39  2001-01-02 05:34:42
## 42  2001-01-02 06:22:16
## 150 2001-01-04 15:46:53
## 180 2001-01-05 05:30:45
## 203 2001-01-05 17:04:38

ggplot() +
  geom_sf(data = df_cat.sf[df_cat$MAG>4,], size = 0.05) +
  geom_sf(data = nz_cropped, fill=alpha("lightgrey", 0), color = 'green') +
  geom_sf(data = df_cat.sf[df_cat$MAG>6,], size = 0.5, color='orange') +
  geom_sf(data = df_cat.sf[df_cat$MAG>7,], size = 0.5, color='red') +
  ggtitle("Map of event locations")+
  geom_rect(aes(xmin=longLims[1], xmax=longLims[2], ymin=latLims[1], ymax=latLims[2]), color="blue",fill=NA)

ggplot() + 
  geom_point(data=df_cat[df_cat$MAG>4,], aes(dateTime, LATITUDE), size=0.1) +
  geom_point(data=df_cat[df_cat$MAG>6,], aes(dateTime, LATITUDE), size=1.2, color='orange') +
  geom_point(data=df_cat[df_cat$MAG>7,], aes(dateTime, LATITUDE), size=1.5, color='red') +
  ggtitle("New Zealand latitude-time plot with Northern Island highlighted") +
  geom_rect(aes(xmin = as.POSIXct(ymd("2001-1-1")), xmax = as.POSIXct(ymd("2012-1-1")), ymin = latLims[1], ymax = latLims[2]), alpha = 0.4, fill='blue')

ggplot() + 
  geom_point(data=df_cat.north[df_cat.north$MAG>4,], aes(dateTime, LATITUDE), size=0.1) +
  geom_point(data=df_cat.north[df_cat.north$MAG>6,], aes(dateTime, LATITUDE), size=1.2, color='orange') +
  geom_point(data=df_cat.north[df_cat.north$MAG>7,], aes(dateTime, LATITUDE), size=1.5, color='red') +
  ggtitle("New Zealand latitude-time plot with Northern Island highlighted") +
  geom_rect(aes(xmin = as.POSIXct(ymd("2001-1-1")), xmax = as.POSIXct(ymd("2012-1-1")), ymin = latLims[1], ymax = latLims[2]), alpha = 0.4, fill='blue')

ggplot(df_cat.north, aes(x=dateTime, y=MAG)) + 
  geom_point(size = 0.1) +
  ggtitle("Northern Island magnitude timeseries for M>4") +
  geom_rect( aes(xmin = as.POSIXct(startDate), xmax = as.POSIXct(endDate), ymin = minMAG, ymax = max(df_cat.north$MAG+0.2)), alpha = 0.4, fill=NA, color="blue" )

## Warning: Use of `df_cat.north$MAG` is discouraged. Use `MAG` instead.

minMag <- 4
maxMag <- max(df_cat.north$MAG)

mags <- df_cat.north[df_cat.north$MAG>=minMag,]$MAG

tmp <- hist(mags, breaks=seq(minMag-0.05,maxMag+0.1,0.1), plot=FALSE)

N.counts <- length( tmp$counts)
tmp$cumulativeCounts <- cumsum(tmp$counts[N.counts:1])[N.counts:1]

m.min <- 4
bin_m.min <- which(tmp$mids==m.min)
freq_m.min <- tmp$counts[bin_m.min]
b <- 1.1
x <- tmp$mids
y <- freq_m.min * 10^(-b*(x-m.min))
y.cum <- tmp$cumulativeCounts[bin_m.min] * 10^(-b*(x-m.min))

ggplot() +
  geom_point( aes(x=tmp$mids, y=tmp$counts) ) +
  geom_point( aes(x=tmp$mids, y=tmp$cumulativeCounts) , color='red', pch="+") +
  scale_y_log10() +
  ggtitle(paste("Frequency-magnitude plot with arbitary GR dist: b =", b)) +
  xlab("Magnitude") +
  ylab("log10(Frequency)") +
  geom_line(aes(x=x, y=y)) +
  geom_line(aes(x=x, y=y.cum), color='red') +
  geom_vline( xintercept=m.min, lty=2 )

## Warning: Transformation introduced infinite values in continuous y-axis

b.stability.list <- c()
b.error.list <- c()
m.mean <- c()

max.index.x <- length(x)-5

for( i in 1:max.index.x ){
  mag.threshold <- x[i]
  m.mean[i] <- mean( mags[mags > mag.threshold], na.rm=TRUE )
  b.stability.list[i] <- log10( exp(1) ) / (m.mean[i] - mag.threshold - 0.05)
  b.error.list[i] <- 2.3 * b.stability.list[i]^2 * sd(mags[mags > mag.threshold]) / (sqrt(N.events*(N.events-1)))
}

ggplot() +
  geom_line( aes(x=x[1:max.index.x], y=b.stability.list) ) +
  geom_line( aes(x=x[1:max.index.x], y=b.stability.list+b.error.list), color=2, lty=2 ) +
  geom_line( aes(x=x[1:max.index.x], y=b.stability.list-b.error.list), color=2, lty=2 ) +
  xlab("Magnitude threshold") +
  ylab("Aki b-value estimate") +
  geom_hline(yintercept = 1, lty=3) +
  ggtitle("b-value stability plot for Northern Island")

Extract the 2009 Fiordland Sequence

UTC time 2009-07-15 09:22:29 ISC event 15157724 USGS-ANSS ComCat Local date 15 July 2009 Local time 9:22 pm (NZST) Magnitude 7.8 Mw Depth 12.0 km (7.5 mi) Epicentre 45.762°S 166.562°ECoordinates: 45.762°S 166.562°E Areas affected New Zealand Max. intensity VI (Strong) Tsunami 2.3 m (7 ft 7 in) Aftershocks >100 Casualties 0

eventDate <- ymd_hms( '2009-07-15 09:22:29' )
endDate <- eventDate + days(100)
startDate <- eventDate - days(50)
deltaLat <- 2.4
latLims <- c( -45.762-deltaLat, -45.762+deltaLat)
longLims <- c( 166.562-deltaLat, 166.562+deltaLat)

minMAG <- 4

# Subset the main catalogue
df_cat.subset <- df_cat[df_cat$MAG >= minMAG, ]
df_cat.subset <- df_cat.subset[ (df_cat.subset$LATITUDE>latLims[1]), ]
df_cat.subset <- df_cat.subset[ (df_cat.subset$LATITUDE<latLims[2]), ]
df_cat.subset <- df_cat.subset[ (df_cat.subset$LONGITUDE>longLims[1]), ]
df_cat.subset <- df_cat.subset[ (df_cat.subset$LONGITUDE<longLims[2]), ]

head(df_cat.subset)

##        YR MO DY HR MN   SEC LATITUDE LONGITUDE  DEPTH MAG NO RMSRES      x
## 163  2001  1  4 20  9 58.04 -45.0396  167.3916  91.06 4.0 18   0.14 -88.59
## 1067 2001  1 22 20 41 24.25 -44.9499  167.4568 101.89 4.6 21   0.09 -91.82
## 1552 2001  2  1  5  2 29.32 -44.5924  168.1326  76.88 4.0 21   0.13 -80.52
## 2190 2001  2 12  3  6 43.91 -46.7453  165.1765  12.75 4.3  7   0.35 -76.92
## 2194 2001  2 12  4 34 17.26 -46.8243  165.3218  30.00 5.0 11   0.38 -62.82
## 3122 2001  2 25 20  1 28.97 -44.6792  168.1733  11.52 4.3 23   0.14 -71.54
##           y GAP DMIN RZDM NP NS SE_OT  SE_H  SE_Z Q event_num       date
## 163  569.18 187   52  1.8 12  6  0.08  1.00  2.00 C       163 2001-01-04
## 1067 558.42 191   65  1.6 13  8  0.06  1.00  2.00 C      1067 2001-01-22
## 1552 492.64 209   86  0.9 14  7  0.05  1.00  2.00 C      1552 2001-02-01
## 2190 825.47 315  244  0.1  6  1  1.03 21.36 26.15 E      2190 2001-02-12
## 2194 823.67 319  215  0.1  7  4  0.71  8.41 16.22 E      2194 2001-02-12
## 3122 497.43 129   20  0.6 17  6  0.02  0.25  0.64 A      3122 2001-02-25
##                 dateTime
## 163  2001-01-04 20:09:58
## 1067 2001-01-22 20:41:24
## 1552 2001-02-01 05:02:29
## 2190 2001-02-12 03:06:43
## 2194 2001-02-12 04:34:17
## 3122 2001-02-25 20:01:28

ggplot() +
  geom_sf(data = df_cat.sf[df_cat$MAG>4,], size = 0.05) +
  geom_sf(data = nz_cropped, fill=alpha("lightgrey", 0), color = 'green') +
  geom_sf(data = df_cat.sf[df_cat$MAG>6,], size = 0.5, color='orange') +
  geom_sf(data = df_cat.sf[df_cat$MAG>7,], size = 0.5, color='red') +
  ggtitle("Map of event locations")+
  geom_rect(aes(xmin=longLims[1], xmax=longLims[2], ymin=latLims[1], ymax=latLims[2]), color="blue",fill=NA)

ggplot() + 
  geom_point(data=df_cat[df_cat$MAG>4,], aes(dateTime, LATITUDE), size=0.1) +
  geom_point(data=df_cat[df_cat$MAG>6,], aes(dateTime, LATITUDE), size=1.2, color='orange') +
  geom_point(data=df_cat[df_cat$MAG>7,], aes(dateTime, LATITUDE), size=1.5, color='red') +
  ggtitle("New Zealand latitude-time plot") +
  geom_rect(aes(xmin = as.POSIXct(ymd("2001-1-1")), xmax = as.POSIXct(ymd("2012-1-1")), ymin = latLims[1], ymax = latLims[2]), alpha = 0.4, fill='blue')

ggplot() + 
  geom_point(data=df_cat[df_cat$MAG>4,], aes(dateTime, LATITUDE), size=0.1) +
  geom_point(data=df_cat[df_cat$MAG>6,], aes(dateTime, LATITUDE), size=1.2, color='orange') +
  geom_point(data=df_cat[df_cat$MAG>7,], aes(dateTime, LATITUDE), size=1.5, color='red') +
  ggtitle("New Zealand latitude-time plot") +
  geom_rect(aes(xmin = as.POSIXct(startDate), xmax = as.POSIXct(endDate), ymin = latLims[1], ymax = latLims[2]), alpha = 0.4, fill='blue')

ggplot(df_cat.subset, aes(x=dateTime, y=MAG)) + 
  geom_point(size = 0.1) +
  ggtitle("New Zealand magnitude timeseries for M>4") +
  geom_rect( aes(xmin = as.POSIXct(startDate), xmax = as.POSIXct(endDate), ymin = minMAG, ymax = max(df_cat.subset$MAG+0.2)), alpha = 0.4, fill=NA, color="blue" )

## Warning: Use of `df_cat.subset$MAG` is discouraged. Use `MAG` instead.

minMag <- 4
maxMag <- max(df_cat.subset$MAG)

mags <- df_cat.subset[df_cat.subset$MAG>=minMag,]$MAG

tmp <- hist(mags, breaks=seq(minMag-0.05,maxMag+0.1,0.1), plot=FALSE)

N.counts <- length( tmp$counts)
tmp$cumulativeCounts <- cumsum(tmp$counts[N.counts:1])[N.counts:1]

m.min <- 4
bin_m.min <- which(tmp$mids==m.min)
freq_m.min <- tmp$counts[bin_m.min]
b <- 1.1
x <- tmp$mids
y <- freq_m.min * 10^(-b*(x-m.min))
y.cum <- tmp$cumulativeCounts[bin_m.min] * 10^(-b*(x-m.min))

ggplot() +
  geom_point( aes(x=tmp$mids, y=tmp$counts) ) +
  geom_point( aes(x=tmp$mids, y=tmp$cumulativeCounts) , color='red', pch="+") +
  scale_y_log10() +
  ggtitle(paste("Frequency-magnitude plot with arbitary GR dist: b =", b)) +
  xlab("Magnitude") +
  ylab("log10(Frequency)") +
  geom_line(aes(x=x, y=y)) +
  geom_line(aes(x=x, y=y.cum), color='red') +
  geom_vline( xintercept=m.min, lty=2 )

## Warning: Transformation introduced infinite values in continuous y-axis

b.stability.list <- c()
b.error.list <- c()
m.mean <- c()

max.index.x <- length(x)-5

for( i in 1:max.index.x ){
  mag.threshold <- x[i]
  m.mean[i] <- mean( mags[mags > mag.threshold], na.rm=TRUE )
  b.stability.list[i] <- log10( exp(1) ) / (m.mean[i] - mag.threshold - 0.05)
  b.error.list[i] <- 2.3 * b.stability.list[i]^2 * sd(mags[mags > mag.threshold]) / (sqrt(N.events*(N.events-1)))
}

ggplot() +
  geom_line( aes(x=x[1:max.index.x], y=b.stability.list) ) +
  geom_line( aes(x=x[1:max.index.x], y=b.stability.list+b.error.list), color=2, lty=2 ) +
  geom_line( aes(x=x[1:max.index.x], y=b.stability.list-b.error.list), color=2, lty=2 ) +
  xlab("Magnitude threshold") +
  ylab("Aki b-value estimate") +
  geom_hline(yintercept = 1, lty=3) +
  ggtitle("b-value stability plot for catalogue subset")